
setwd("c:/users/james/desktop/dropbox/transparency_and_democracy/transparencyindex/PAReplicationMaterials/regressions/")

library(R2jags)
library(foreign)

load("TransparencyABPrep.RData")

attach.jags(results)

HRV.sigma.burqual<-c(1:2000)
HRV.sigma.corrup<-c(1:2000)
HRV.sigma.laword<-c(1:2000)

HRV.sigma.burqual<- sigma.burqual
HRV.sigma.corrup <- sigma.corrup
HRV.sigma.laword <- sigma.laword

HRV.fit<- data.frame(HRV.sigma.burqual, HRV.sigma.corrup, HRV.sigma.laword)

write.dta(HRV.fit, file="HRVModelFit.dta")

rm(list=ls())

load("FracABPResults.RData")

attach.jags(results.frac)

HRV.fit<-read.dta("HRVModelFit.dta")

corrupFrac<-hist(sigma.corrup)
corrupHRV<- hist(HRV.fit$HRV_sigma_corrup)

burqualFrac<-hist(sigma.burqual)
burqualHRV<-hist(HRV.fit$HRV_sigma_burqual)

lawordFrac<-hist(sigma.laword)
lawordHRV<-hist(HRV.fit$HRV_sigma_laword)

pdf('CorruptionModelFit.pdf')
plot(corrupFrac, col=rgb(0,0,1,1/4), freq=FALSE, xlab="sigma", main="Corruption Model Error", ylim=c(0,15))
plot(corrupHRV, col=rgb(1,0,1,1/4), freq=FALSE, add=TRUE)
dev.off()

pdf('BurQualModelFit.pdf')
plot(burqualFrac, col=rgb(0,0,1,1/4), freq=FALSE, xlab="sigma", main="Bureaucratic Quality Model Error", ylim=c(0,20))
plot(burqualHRV, col=rgb(1,0,1,1/4), freq=FALSE, add=TRUE)
dev.off()

pdf('LawOrdModelFit.pdf')
plot(lawordFrac, col=rgb(0,0,1,1/4), freq=FALSE, xlab="sigma", main="Law and Order Model Error")
plot(lawordHRV, col=rgb(1,0,1,1/4), freq=FALSE, add=TRUE)
dev.off()

fracmargeffectCorrup.sim <- array(NA, c(2000,100))
fracmargeffectLawOrd.sim <- array(NA, c(2000,100))
fracmargeffectBurQual.sim <- array(NA, c(2000,100))

polity.sim<- -10

for (i in 1:100){
  fracmargeffectCorrup.sim[, i]  <-  gamma.corrup[,4]*sd(frac_reported) + gamma.corrup[,5]*sd(frac_reported)*polity.sim
  fracmargeffectLawOrd.sim[, i]  <-  gamma.laword[,4]*sd(frac_reported) + gamma.laword[,5]*sd(frac_reported)*polity.sim
  fracmargeffectBurQual.sim[, i]  <-  gamma.burqual[,4]*sd(frac_reported) + gamma.burqual[,5]*sd(frac_reported)*polity.sim
  polity.sim<- polity.sim + .2
  }

fraccorrupEst <- c(1:100)

fraccorrupLB <- c(1:100)
fraccorrupUB <- c(1:100)

fraccorrupHPD <- cbind(fraccorrupLB, fraccorrupUB)

fraclawordEst <- c(1:100)

fraclawordLB <- c(1:100)
fraclawordUB <- c(1:100)

fraclawordHPD <- cbind(fraclawordLB, fraclawordUB)

fracburqualEst <- c(1:100)

fracburqualLB <- c(1:100)
fracburqualUB <- c(1:100)

fracburqualHPD <- cbind(fracburqualLB, fracburqualUB)

for (i in 1:100){
  fraccorrupEst[i] <- mean(fracmargeffectCorrup.sim[,i])
  fraclawordEst[i] <- mean(fracmargeffectLawOrd.sim[,i])
  fracburqualEst[i] <- mean(fracmargeffectBurQual.sim[,i])

  fraccorrupHPD[i,] <- HPDinterval(as.mcmc(as.matrix(fracmargeffectCorrup.sim[,i])))
  fraclawordHPD[i,] <- HPDinterval(as.mcmc(as.matrix(fracmargeffectLawOrd.sim[,i])))
  fracburqualHPD[i,] <- HPDinterval(as.mcmc(as.matrix(fracmargeffectBurQual.sim[,i])))
  }

polity.sim<-c(1:100)
polity.sim<-(polity.sim-51)/5

pdf('FracCorruptEffect.pdf')
plot(x=c(), y=c(), xlim=c(-10, 10), ylim=c(-1,1), xlab="Polity2", ylab="Marginal Corruption Effect", main="Fraction Reported and Corruption")
lines(y=fraccorrupEst, x=polity.sim)
lines(y=fraccorrupHPD[,1], x=polity.sim, lty=2)
lines(y=fraccorrupHPD[,2], x=polity.sim, lty=2)
abline(h=0, lty=3)
dev.off()


pdf('FraclawordEffect.pdf')
plot(x=c(), y=c(), xlim=c(-10, 10), ylim=c(-1,1), xlab="Polity2", ylab="Marginal Corruption Effect", main="Fraction Reported and Law and Order")
lines(y=fraclawordEst, x=polity.sim)
lines(y=fraclawordHPD[,1], x=polity.sim, lty=2)
lines(y=fraclawordHPD[,2], x=polity.sim, lty=2)
abline(h=0, lty=3)
dev.off()

pdf('FracburqualEffect.pdf')
plot(x=c(), y=c(), xlim=c(-10, 10), ylim=c(-1,1), xlab="Polity2", ylab="Marginal Corruption Effect", main="Fraction Reported and Bureaucratic Quality")
lines(y=fracburqualEst, x=polity.sim)
lines(y=fracburqualHPD[,1], x=polity.sim, lty=2)
lines(y=fracburqualHPD[,2], x=polity.sim, lty=2)
abline(h=0, lty=3)
dev.off()

write.dta(data.frame(fracmargeffectCorrup.sim), file="FracMargEffectCorrupt.dta")
write.dta(data.frame(fracmargeffectBurQual.sim), file="FracMargEffectBurQual.dta")
write.dta(data.frame(fracmargeffectLawOrd.sim), file="FracMargEffectLawOrd.dta")


rm(list=ls())

load("TransparencyABPRep.RData")

FracMargEffectCorrup<-read.dta("FracMargEffectCorrupt.dta")
FracMargEffectBurQual<-read.dta("FracMargEffectBurQual.dta")
FracMargEffectLawOrd<-read.dta("FracMargEffectLawOrd.dta")

attach.jags(results)

k<-0
transparency.means<-c(1:2120)

for (i in 1:106){
  for (j in 1:20){
  k<-k+1
  transparency.means[k] <- mean(transparency[,i,j])
  }
}

HRVmargeffectCorrup.sim <- array(NA, c(2000,100))
PapermargeffectCorrup.sim<- array(NA, c(2000,100))
HRVmargeffectLawOrd.sim <- array(NA, c(2000,100))
PapermargeffectLawOrd.sim<- array(NA, c(2000,100))
HRVmargeffectBurQual.sim <- array(NA, c(2000,100))
PapermargeffectBurQual.sim<- array(NA, c(2000,100))

polity.sim<- -10

for (i in 1:100){
  HRVmargeffectCorrup.sim[, i]  <-  gamma.corrup[,4]*sd(transparency.means) + gamma.corrup[,5]*sd(transparency.means)*polity.sim
  PapermargeffectCorrup.sim[,i] <- gamma.corrup[,6]*sd(newscirc) + gamma.corrup[,7]*sd(newscirc)*polity.sim
  HRVmargeffectLawOrd.sim[, i]  <-  gamma.laword[,4]*sd(transparency.means) + gamma.laword[,5]*sd(transparency.means)*polity.sim
  PapermargeffectLawOrd.sim[,i] <- gamma.laword[,6]*sd(newscirc) + gamma.laword[,7]*sd(newscirc)*polity.sim
  HRVmargeffectBurQual.sim[, i]  <-  gamma.burqual[,4]*sd(transparency.means) + gamma.burqual[,5]*sd(transparency.means)*polity.sim
  PapermargeffectBurQual.sim[,i] <- gamma.burqual[,6]*sd(newscirc) + gamma.burqual[,7]*sd(newscirc)*polity.sim
  polity.sim<- polity.sim + .2
  }

DiffMargEffectCorrup<-array(NA, c(2000,100))
DiffMargEffectBurQual<-array(NA, c(2000,100))
DiffMargEffectLawOrd<-array(NA, c(2000,100))

DiffMargEffectCorrup <- HRVmargeffectCorrup.sim - FracMargEffectCorrup
DiffMargEffectBurqual <- HRVmargeffectBurQual.sim - FracMargEffectBurQual
DiffMargEffectLawOrd <- HRVmargeffectLawOrd.sim - FracMargEffectLawOrd

diffcorrupEst <- c(1:100)

diffcorrupLB <- c(1:100)
diffcorrupUB <- c(1:100)

diffcorrupHPD <- cbind(diffcorrupLB, diffcorrupUB)

difflawordEst <- c(1:100)

difflawordLB <- c(1:100)
difflawordUB <- c(1:100)

difflawordHPD <- cbind(difflawordLB, difflawordUB)

diffburqualEst <- c(1:100)

diffburqualLB <- c(1:100)
diffburqualUB <- c(1:100)

diffburqualHPD <- cbind(diffburqualLB, diffburqualUB)

for (i in 1:100){
  diffcorrupEst[i] <- mean(DiffMargEffectCorrup[,i])
  difflawordEst[i] <- mean(DiffMargEffectLawOrd[,i])
  diffburqualEst[i] <- mean(DiffMargEffectBurqual[,i])

  diffcorrupHPD[i,] <- HPDinterval(as.mcmc(as.matrix(DiffMargEffectCorrup[,i])))
  difflawordHPD[i,] <- HPDinterval(as.mcmc(as.matrix(DiffMargEffectLawOrd[,i])))
  diffburqualHPD[i,] <- HPDinterval(as.mcmc(as.matrix(DiffMargEffectBurqual[,i])))
  }

polity.sim<-c(1:100)
polity.sim<-(polity.sim-51)/5

pdf('DiffCorruptEffect.pdf')
plot(x=c(), y=c(), xlim=c(-10, 10), ylim=c(-1,1), xlab="Polity2", ylab="HRV Effect - Frac Reported Effect", main="Difference in Marginal Effects on Corruption")
lines(y=diffcorrupEst, x=polity.sim)
lines(y=diffcorrupHPD[,1], x=polity.sim, lty=2)
lines(y=diffcorrupHPD[,2], x=polity.sim, lty=2)
abline(h=0, lty=3)
dev.off()


pdf('DifflawordEffect.pdf')
plot(x=c(), y=c(), xlim=c(-10, 10), ylim=c(-1,1), xlab="Polity2", ylab="HRV Effect - Frac Reported Effect", main="Difference in Marginal Effects on Law and Order")
lines(y=difflawordEst, x=polity.sim)
lines(y=difflawordHPD[,1], x=polity.sim, lty=2)
lines(y=difflawordHPD[,2], x=polity.sim, lty=2)
abline(h=0, lty=3)
dev.off()

pdf('DiffburqualEffect.pdf')
plot(x=c(), y=c(), xlim=c(-10, 10), ylim=c(-1,1), xlab="Polity2", ylab="HRV Effect - Frac Reported Effect", main="Difference in Marginal Effects on Bureaucratic Quality")
lines(y=diffburqualEst, x=polity.sim)
lines(y=diffburqualHPD[,1], x=polity.sim, lty=2)
lines(y=diffburqualHPD[,2], x=polity.sim, lty=2)
abline(h=0, lty=3)
dev.off()
